function [sresults] = estimate_rcnl_supply( ...
    df,          ...
    cdindex,     ...
    dresults,    ...
    derMat,      ...
    IVDtilda,    ...
    instruments  ...
)


    [~, w, ~, ~, fess, ~, IVS_noFE] = final_variables(df, instruments, 'og');

    options = optimset(         ...
        'GradObj',     'off',   ...
        'MaxIter',     500000,  ...
        'MaxFunEvals', 1000000, ...
        'Display',     'iter',  ...
        'TolFun',      1e-4,    ...
        'TolX',        1e-4     ...
    );

    IVStilda   = IVS_noFE - fess * ((fess' * fess) \ (fess' * IVS_noFE));
    invAStilda = inv(IVStilda' * IVStilda);

    %%%%%%%%%%%%%%%%
    %  First Step  %
    %%%%%%%%%%%%%%%%

    tic
    f_obj = @(sparams) rcnl_gmmobj_supply( ...
        sparams,    ...
        cdindex,    ...
        df.cdid,    ...
        df.firmid,  ...
        df.prodid,  ...
        df.yearid,  ...
        df.bmc,     ...
        df.s_jt,    ...
        df.p_jt,    ...
        w,          ...
        derMat,     ...
        IVStilda,   ...
        invAStilda, ...
        fess        ...
    );

    [sparams, sfval, sexit, ~] = fminsearch(f_obj, [0.5], options);
    [sfval, gamma, eta] = f_obj(sparams);
    shours = toc / 60 / 60;
    disp(['Elapsed: ' num2str(shours)]);

    sresults         = struct;
    sresults.eta     = eta;
    sresults.sfval   = sfval;
    sresults.gamma   = gamma;
    sresults.kappa   = sparams;
    sresults.ngamma  = size(w, 2) - size(fess, 2);
    sresults.sparams = [sparams; sresults.gamma(1:sresults.ngamma)];


end
